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Abstract 

A fast multigrid solver for the steady incom- 
pressible Euler equations is presented. Unlike time- 
marching schemes, this approach uses relaxation of 
the steady equations. Application of this method 
results in a discretization that correctly distin- 
guishes between the a dvection and elliptic parts 
of the operator, allowing efficient smoothers to be 
constructed. Solvers for both unstructured trian- 
gular grids and structured quadrilateral grids have 
been written. Computations for channel Bow and 
Bow over a nonlifting airfoil have computed. Us- 
ing Gauss-Seidel relaxation ordered in the flow di- 
rection, textbook multigrid convergence rates of 
nearly one order-of-magnitude residual reduction 
per multigrid cycle are achieved, independent of the 
grid spacing. This approach also may be applied to 
the compressible Euler equations and the incom- 
pressible Navier-Stokes equations. 

Introduction 

One of the critical needs in computational fluid dy- 
namics is faster flow solvers. Multigrid is a well known 
method of convergence acceleration which is widely used 
in Euler and Reynolds-averaged Navier-Stokes codes. 
These applications of multigrid generally are based on 
the unsteady equations using some temporal integrator 
as the smoother, combined with a full-approximation 
scheme (FAS) multigrid iteration. A common approach 
is one originally proposed by Jameson^. Starting with 
the unsteady equations, a finite-volume spatial dis- 
cretization with explicit artificial viscosity is combined 
with a Runge-Kutta time integration as a smoother. An 
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alternative approach ’ ’ is to use upwind-differencing 

and implicit time integration as the smoother. However, 

Copyright ©1997 by the American Institute of Aeronau- 
tics and Astronautics, Inc. No copyright is asserted in the 
United States under Title 17, U.S. Code. The U.S. Govern- 
ment has a royalty-free license to exercise all rights under the 
copyright claimed herein for government purposes. All other 
rights are reserved by the copyright owner. 

* Research Scientist 

1 Senior Staff Scientist 


these approaches have resulted in poor multigrid effi- 
ciency. When applied to high Reynolds number flows 
over complex geometries, convergence rates are often 
worse than 0.99. There is clearly a need to develop sub- 
stantially more efficient multigrid solvers. 
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According to Brandt , one of the major obstacles 
to achieving better multigrid performance for advection 
dominated flows is that the coarse grid provides only a 
fraction of the needed correction for smooth error com- 
ponents. This obstacle can be removed by designing a 
solver that effectively distinguishes between the elliptic, 
parabolic, and hyperbolic (advection) factors of the sys- 
tem and treats each one appropriately. For instance, 
advection can be treated by space marching, while el- 
liptic factors can be treated by multigrid. The effi- 
ciency of such an algorithm will be essentially identi- 
cal to that of the solver for the elliptic factor only, and 
thereby attain so-called “textbook” multigrid efficiency. 
Brandt presents an approach called “distributive relax- 
ation” by which one can construct smoothers that ef- 
fectively distinguish between the different factors of the 
operator. Using this approach, Brandt and Yavneh have 
demonstrated textbook multigrid for the incompressible 
Navier-Stokes equations 0 . Their results are for a simple 
geometry and a Cartesian grid, using a staggered-grid 
discretization of the equations. 

In a closely related approach, Ta’asan presents a fast 
multigrid solver for the compressible Euler equations. 
This method is based on a set of “canonical variables” 
which express the steady Euler equations in terms of an 

O 

elliptic and a hyperbolic partition 0 . Ta’asan uses this 
partition to guide the discretization of the equations. A 
staggered grid is used, with different variables residing at 
cell, vertex, and edge centers. In Reference 7 it is shown 
that ideal multigrid efficiency can be achieved for the 
compressible Euler equations for two-dimensional sub- 
sonic flow using body-fitted grids. One possible limita- 
tion of the use of canonical variables is that the partition 
of the inviscid equations is not directly applicable to the 
viscous equations. 

In this paper, an alternative to distributive relax- 
ation and to Ta’asan’s canonical variable decomposi- 
tion is presented. It is a generalization of the approach 
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of Sidilkover and Ascher^, and is described in detail 

by Sidilkover^ . This approach can be classified as a 
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method of the Weighted Gauss-Seidel type . A con- 
ventional vertex-based finite-volume or finite-difference 
discretization of the primitive variables is used, avoid- 
ing the need for staggered grids. This simplifies the re- 
striction and prolongation operations, because the same 
operator can be used for all variables. A projection op- 
erator is applied to the system of equations, resulting 
in a Poisson equation for the pressure. By applying the 
projection operator to the discrete equations rather than 
to the differential equations, the proper boundary con- 
dition on the pressure is satisfied directly. The Poisson 
equation for the pressure may be treated by Gauss-Seidel 
relaxation, while the advection terms of the momentum 
equation are treated by space-marching. Because the el- 
liptic and advection parts of the system are decoupled, 
ideal multigrid efficiency can be achieved. Compared to 
distributive relaxation and the canonical variables ap- 
proaches, this method is extremely simple. 

Formulation of the Problem 

The incompressible Euler equations in primitive vari- 
ables are 

UU X + VUy + p x = 0 , 

UV X + Wy + Py = 0, 

U x + Vy = 0, 

where u and v are the components of the velocity in 
the x and y directions, respectively, andp is the pressure. 
The density is taken to be one. The advection operator 
is defined by 

Q = ud x + vdy , (1) 

where d x , d y are the partial differentiation operators. 
The Euler equations may be written as 

/ Q 0 d x \ / u\ 

Lq = 0 Q d y r = 0. (2) 

\d x dy 0 J \pj 

Introducing the adjoint to Q, defined by 

<?*(/) = —d x (uf) - dy(vf), (3) 

a projection operator P is defined: 

// 0 0\ 

P = 0 / 0 . (4) 

\d x d y Q*J 

Applying the projection operator to the Euler equations 
yields 

/ Q 0 d x \ /u\ 

Lq = PLq =10 Q d y r + s.p.t., (5) 

\0 0 A / \pj 

where A is the Laplacian. The matrix operator on 
the right-hand side consists of the principal part of L, 


and “s.p.t.” are the subprincipal terms, in the termi- 
nology of Brandt . These terms arise because the co- 
efficients u and v in the operators Q and Q * are not 
constant. It is important to note that the subprincipal 
terms can be ignored for the purpose of constructing a 
relaxation scheme. 

A feature of Eq. (5) is that the correct boundary con- 
dition for the pressure at a solid wall is an automatic 
result of the tangency boundary condition on the veloc- 
ity. The equation for the pressure at the wall reduces 
to 

( — udyV + vdyii + d x p)n x 

+ ( ud x v — vd x u + d y p)n y = 0, (6) 

where n x , h y are the components of the unit normal 
at the wall. Because the normal component of veloc- 
ity, uh x + vh y , is zero, it can be shown that Eq. (6) 
reduces to the momentum equation in the normal direc- 
tion. No auxiliary boundary condition on the pressure 
is needed. 

The operator on the left-hand side of Eq. (5) is up- 
per triangular. The pressure satisfies a Poisson equa- 
tion for which a conventional relaxation method, such 
as Gauss-Seidel, can be applied. Upwind differencing of 
the advection operator in the momentum equations al- 
lows downstream relaxation to be used. The strategy 
used to relax the system is to first update the pressure. 
The pressure update contributes to the velocity update 
through the gradient terms in the right-hand column of 
the operator in Eq. (5). Finally, the velocity components 
are updated by relaxing in the stream direction. 

Solution Procedure 

Rather than discretizing Eq. (5) directly, the first 
step of numerical solution procedure is to discretize 
Eq. (2). The relaxation scheme is constructed by 
applying the projection operator P at the discrete 
level rather than the differential level. A sequence of 
grids Gk, Gk— i, • • • , Go is used, where Gk is the finest 
and Go the coarsest. Let La, be the discrete approxima- 
tion to the operator L and qA be the solution on the fc-th 
grid. This system has the form LAqA = fi, where the en- 
tries of La are 3x3 block matrices which operate on the 
unknowns ( u,v,p) T at each grid vertex. A general iter- 
ation scheme is constructed by writing the operator La 
as La = Ma — Na, where the splitting is chosen such 
that Ma is easily inverted. Lexicographic Gauss-Seidel 
is obtained by taking Ma to be the block lower-triangular 
matrix resulting from ignoring the blocks above the di- 
agonal of La. A further simplification is obtained if the 
diagonal blocks of Ma contain only those entries corre- 
sponding to the principal part of the operator. Because 
the operator in Eq. (5) is upper triangular, the diagonal 
blocks will then be 3 x 3 upper triangular matrices. 

Letting q£ be the n- th iterate of the solution, the 
iteration is 

MAqA + 1 =fA + NAq£. 
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Because the operator L& is nonlinear, M& and N& will 
be functions of qjl and q£ +1 . Letting 5qjl = q)! +1 — qjl, 
the iteration may be rewritten as 

M .'q. f. I. q. • (7) 

Because M& is block lower-triangular, S qjl is found by 
forward substitution and inverting a 3 x 3 diagonal block 
at each vertex. The diagonal blocks are upper triangular 
and are easily inverted. 

If upwind differences are used for the advection op- 
erator ( 1 ) and the grid points are ordered in the flow 
direction, then the 3x3 blocks of N& will have zeroes 
in the first two rows. In this case, lexicographic Gauss- 
Seidel relaxation is equivalent to space-marching of the 
advection terms. The advected error is effectively elimi- 
nated in one relaxation sweep and the convergence rate 
becomes that of the Poisson equation for the pressure. It 
is possible to get ideal multigrid convergence rates for the 
system because each component of the error is treated 
appropriately. 

A straightforward FAS multigrid iteration is applied 
to the system of equations. Let Lfc_i be the coarse grid 
operator, I^_ 1 be the hne-to-coarse grid restriction op- 
erator, and /£ -1 be the coarse-to-hne grid prolongation 
operator. If is the current solution on grid k, the 
residual on this grid is r* = fk — L%q%. This leads to the 
coarse-grid equation 

Lfc-iqfc-i = ffc-i = /fc-ii’fc + Lfc-i (fl-iqfc) . (8) 

After solving the coarse-grid equation for q/,--i, the fine- 
grid solution is corrected by 

q* <- q k + /fc _1 (qfc-i - Ik-Ak^j ■ (9) 

■liquation (8) is solved by applying the same relaxation 
procedure that, is used to solve the fine-grid equation. 
Multigrid is applied recursively to the coarse-grid equa- 
tion. On the coarsest grid, many relaxation sweeps are 
performed to insure that, the equation is solved com- 
pletely. A conventional V-cycle or IT-cycle is used. 

Results 

Two flow solvers have been written to implement the 
solution algorithm presented in the previous section. 
One is an unstructured grid code based on triangular 
elements. The program is a modification of the hnit.e- 
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volume code FUN2D ’ “. The discretization of the mo- 

mentum equations uses the advection scheme proposed 
13 

by Giles, $t al. The residuals of the momentum equa- 
tions are computed by a trapezoidal rule integration 
around each triangle of the grid. These residuals are dis- 
tributed to the nodes of the triangle using upwind-biased 
weights that, depend on the local flow direction. One can 
think of this as a rot.at.ed-difference stencil. The Poisson 
equation for the pressure is discretized using a standard 
central-difference approximation for the Laplacian. 

The other code is a structured, quadrilateral grid 
solver. A cell-vertex, finite-volume method is used at 
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Figure 1. Stencil for structured quadrilateral grid. 
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Figure 2. Channel geometry. 


each grid point. The hnit.e-volume cell associated with 
each vertex of the grid is the corresponding dual grid cell, 
as shown in Fig. 1. For vertex E the dual cell is ABCD. 
The advection terms are discretized by a conventional 
one-sided,, second-order difference approximation. The 
gradients of u, v, and p appearing in the pressure Pois- 
son equation are computed by applying Green’s theorem. 
The gradient on the cell face CD of Fig. 1 is computed 
by integrating around the path ECFD. Derivatives of 
these gradients are also evaluated by Green’s theorem. 
This approach yields the standard five-point approxima- 
tion to the Laplacian operator on a uniform grid. 

Solutions for incompressible, inviscid how in a channel 
have been obtained with both solvers. The channel ge- 
ometry and boundary conditions are shown in Fig. 2. 
The shape of the lower wall between 0 < x < 1 is 
y(x) = t sin 2 7T x . For the computations shown here, the 
thickness ratio r is 0.05. The velocity is specihed at the 
inlet and the pressure is specihed at the outlet. At the 
upper and lower walls of the channel, the how t.angency 
condition u-h = 0 is enforced. 

Similar grids were used for both how solvers. Quasi- 
uniform quadrilateral grids were generated. A simple 
shearing transformation was used in the center part of 
the channel to get boundary conforming grids. For the 
unstructured grid solver, the grids were triangulated by 
dividing each quadrilateral cell along a diagonal. A series 
of nested coarse grids was obtained by coarsening the hne 
grids by a factor of two in each coordinate direction. In 
all cases shown below, the coarsest grid was 7x3 vertices. 
Lexicographic Gauss-Seidel relaxation was used. The or- 
dering of the grid vertices was from the lower-left, to the 
upper-right, of the channel, resulting in downstream re- 
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laxation of the advect.ion operator. A V r (2, 1) multigrid 
cycle was used; that is, two relaxation sweeps were per- 
formed on each grid before restricting to the coarse grid, 
and one relaxation sweep was performed after the coarse- 
grid correction was added to the fine-grid solution. If one 
work unit is taken to be a single Gauss-Seidel sweep on 
the finest grid, a single (2,1) cycle is approximately 
4 work units. 

The computed pressure obtained by the unstructured 
grid flow solver on a grid of 97 x 33 vertices is shown in 
Fig. 3. A total of 5 grid levels was used. The Li norm of 
the residual of each equation for each cycle is shown in 
Fig. 4. The convergence rate is approximately 0.18 per 



Figure 3. Pressure, contour increment Ap = 0.01, 
for an unstructured grid of 97 x 33 vertices. 



Figure 5. Pressure, contour increment Ap = 0.01, 
for a structured grid of 97 x 33 vertices. 



multigrid cycle. For Gauss-Seidel relaxation applied to 
the Laplace equation, one can expect an optimum rate 
of 0.125 for a 1/2, .1 ' multigrid cycle if the restriction 
and prolongation operators are perfect, so the current 
scheme gives very nearly the optimum rate. Also the 
convergence shows no sign of stalling but continues at a 
uniform rate. Changing to a LF-cycle resulted in some 
improvement to the convergence rate per cycle, but re- 
quired more work overall. (In two dimensions, a I'F-cycle 
is about 50% more work than a V-cycle.) 

Very similar results were obtained by the structured 
grid flow solver as shown in Figs. 5 and 6. The conver- 
gence rate of the structured grid code is 0.15 per multi- 
grid cycle, slightly better than that of the unstructured 
grid code and very close to the optimal rate. The better 
performance of the structured grid solver is most likely 
because of better restriction and prolongation operators; 



Figure 6. Convergence history, structured grid 
of 97 x 33 vertices. 

the unstructured flow solver performs bilinear interpola- 
tion using only the locations of a fine-grid vertex and 
the three vertices of the coarse-grid cell containing that 
vertex. 

The present scheme also shows ideal multigrid conver- 
gence independent of the grid spacing. In Fig, 7, the L\ 
norm of the pressure residual is shown for the unstruc- 
tured grid flow solver on a 49 x: 17 fine grid using 4 grid 
levels, up to a 385 x 129 fine grid using 7 grid levels. A 
similar comparison for the structured grid flow solver is 
shown in Fig. 8. This shows that convergence is achieved 
in order n operations. 



Figure 7. Comparison of convergence rates on un- 
structured grids. 


Figure 4. Convergence history, unstructured grid 
of 97 x 33 vertices. 


4 


0 




K(2,l) cycles X^(2,l) cycles 


Figure 8. Comparison of convergence rates on struc- 
tured grids. 



Figure 9. Shear flow, {(-component of velocity, 
contour increment Am = 0.01, unstructured grid 
of 97 x 33 vertices. 


Figure 11. Convergence rate, shear flow in a channel, 
unstructured grid of 97 x 33 vertices. 



Figure 12. Grid generated by perturbing the vertices 
of the 49 x 17 grid. 


The flows shown so far are irrotational. To demon- 
strate that the multigrid performance does not degrade 
if vort.icity is present, a solution for a shear flow in a 
channel was obtained. The unstructured grid flow solver 
was used for this computation. A linear shear was speci- 
fied at the inlet; the velocity was u(y) = (l + y)/2, v = 0. 
The initial condition was a uniform flow in the channel. 
This insured that the vorticity would be introduced into 
the channel only by the relaxation. The {(-velocity com- 
ponent is shown in Fig. 9 and the pressure is shown in 
Fig. 10. The convergence rate shown in Fig, If is seen 
to be essentially the same as in Fig. 4. 

All the calculations shown so far were based on a series 
of nested grids: each coarse-grid vertex corresponds to a 
vertex on the next finer grid. For complex geometries it 
may not practical to generate a series of nested unstruc- 
tured grids, and the performance of the multigrid solver 
may be expected to deteriorate. To show the robust- 


ness of the current method, the triangular grid solver 
was run for a series of non-nested coarse grids. These 
were generated by randomly perturbing the locations of 
the vertices on each of the nested grids independently. 
The perturbed 49 x 17 grid is shown in Fig. 12. The 
computed pressure on a perturbed 97 x, 33 fine grid with 
5 grid levels is shown in Fig. 13 and the convergence 
rate is shown in Fig. 14. The pressure contours are 
very smooth, showing no sign of the lack of grid smooth- 
ness. The asymptotic convergence rate has deteriorated 
to a still-respectable 0.24 per cycle. 

Solutions for incompressible, inviscid, nonlifting flow 
over a Karman-Trefftz airfoil have been obtained with 
the structured grid solver. A cylinder is mapped into a 
Karman-Trefftz airfoil using the conformal mapping 



Figure 10. Shear flow, pressure, contour increment 
Ap = 0.01, unstructured grid of 97 x 33 vertices. 
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Figure 13. Pressure, contour increment Ap = 0.01, 
randomly perturbed unstructured grid of 97 x33 ver- 
tices. 



Figure 14. Convergence rate, randomly perturbed 
unstructured grid of 97 x 33 vertices, 5 grid levels, 
V(2, 1) cycle. 
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Figure 15. Pressure coefficient, nonlifting Karman- 
Trefftz airfoil, 97 x 49 grid. 


where s is the coordinate in the transformed plane, £ is 
the coordinate in the circle plane, and e is the trailing- 
edge angle divided by 7r. A trailing-edge angle of 10° is 
used, resulting in an airfoil of approximately 15% thick- 
ness. The airfoil flow is solved on a finite domain. At 
inflow points along the outer boundary the total pres- 
sure and flow inclination angle are specified. For outflow 
points the pressure is specified. The specified quantities 
are determined from the complex potential function for 
flow past the airfoil. On the airfoil surface the tangency 
condition is enforced. 

A fine grid for the airfoil calculation was constructed 
by generating an O-grid with unit aspect ratio cells in 
the circle plane and mapping it to the airfoil plane. The 
outer boundary of the domain was roughly 13 chords 
from the airfoil. The coarsest grid in the grid sequence 
used for the multigrid solver contains 12 x 6 cells. On 
each grid, relaxation was performed along radial grid 
lines, sweeping from the inflow boundary to the airfoil 
surface over the forward half of the domain, and from the 
airfoil to the outflow boundary over the rearward half 
of the domain. Each sweep started along the stagnation 
streamline, proceeded over the upper half of the domain, 
and then over the lower half of the domain. 

Comparisons between computed and analytic sur- 
face pressure coefficients for nonlifting flow around the 
Karman-Trefftz airfoil are shown in Figs. 15 and 16. 
A W(2, 1) multigrid cycle was used for these computa- 
tions. The computed surface pressure distribution for 
the 97 x 49 grid is in fairly good agreement with the 
analytic one. There are small differences near the lead- 
ing edge and there is the usual difficulty in computing 
full pressure recovery at the trailing edge of the airfoil. 
Note that there is no clustering of the grid in these re- 
gions, which exacerbates the problem. The solution on 
the 193 x 97 grid is in excellent agreement with the 
analytic one, again except for the pressure recovery at 
the trailing edge. The computed pressure contours for 
the 193 x 97 are shown in Fig. 17. 

The convergence history for the 193 x 97 grid is pre- 
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Figure 16. Surfacf. pressure coefficient, nonlifting 
Karman-Trefftz airfoil, 193 x 97 grid. 


sented in Fig. 18. The L\ norm of the residuals for all 
three equations i* shown. The average rate of reduction 
of the residual is 0.12 per W(2, 1). A comparison of the 
convergence rates of the pressure residual for three grid 
densities is shown in Fig. 19. A slight deterioration of 
the convergence rate with increasing grid refinement is 
observed: on the 385 x 197 grid, the rate is 0.15 per cy- 
cles. Nevertheless, as with the channel flow results, the 
convergence rates are very nearly grid independent. 


Extension to Compressible Flow 

The general approach outlined here may be readily ex- 
tended to the compressible Euler equation. In primitive 
variables the equations are 


(Q 

0 

0 

0 ^ 


(s\ 

0 

pQ 

0 

d x 


U 

0 

0 

pQ 
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j?Q) 


W 


( 10 ) 


6 






4r 



.4 t — , — 1 — — 

-6 -4 -2 0 2 4 

X 


Figure 17. Pressure coefficient contours, nonlifting 
Karman-Trefftz airfoil, 193 x 97 grid. 



Figure 18. Convergence history, nonlifting Karman- 
Trefftz airfoil, 193 x 97 vertices. 


where p is the density, c is the speed of sound, and s is 
the entropy. The projection operator P for this system 
is 


(I 0 0 0 \ 

0/00 
0 0/0 
\o d x dy <r) 


CD 


The operators Q and Q * are defined as in Eqs. (1) 
and (3). Applying this to |3q. (10) and ignoring the 
subprincipal terms as before yields 


PLq 


(Q 

0 

0 

0 ^ 


(s\ 

0 

pQ 

0 

d x 


U 

0 

0 

pQ 

By 


V 

V° 

0 

0 

A - M 2 dl y 


W 


+ s.p.t., 

( 12 ) 


where M is the Mach number and d s is the partial deriva- 
tive in the streamwise direction. 

The most significant difference between the compress 1 - 
ible and the incompressible equations is that a Prandtl- 
Glauert-like operator acts on the pressure. Note that 



Figure 19. Comparison of convergence rates for non- 
lifting Karman-Trefftz airfoil. 


this system approaches the system for the incompressible 
equations in the limit of vanishing Mach number. For 
subsonic flow the compressible equations can be solved 
by the same relaxation scheme as the incompressible 
equations. 

Although Eq. (10) is written in primitive variables, 
it is straightforward to use a conservative, finite-volume 
discretization of the system of equations. To relax the 
system, a transformation to primitive variables is made 
to compute the update of the solution values. For 
transonic flows with shocks a stable discretization can 
be achieved using the multidimensional high-resolution 
scheme developed by Sidilkover^ . 

Conclusions 

A multigrid algorithm which yields textbook multigrid 
efficiency for the steady Euler equations has been devel- 
oped. It has the virtue of simplicity; conventional finite- 
difference or finite-volume discretizations of the govern- 
ing equations may be used, allowing flexibility in the 
choice of the underlying numerical method. The correct 
boundary condition for the pressure equation is obtained 
directly. Unlike time-marching approaches* the speed of 
the method does not degrade for low-speed flows, and 
the correct incompressible limit is recovered. Finally, 
this method can be applied to incompressible, viscous 
flow following the ideas of Sidilkover and Ascher^. 
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